C:  main routine to find new potential Phi
C     Subroutine new_phi(r,alpha,theta,dtheta,C,phiu,rho0,
C        			 LogInit,LogRapid);  *
C---- '-->Entering Node %_newphi:'
C-      Subroutine new_phi(r,alpha,theta,dtheta,C,LogInit,LogRapid,rotlaw);

      Subroutine new_phi(C,LogInit,LogRapid,rotlaw)
C-*****************************************************************
      use accur
      use cylinderGrid
      Implicit none
C-        _include accur;
C-        include 'accur.for';
C-       include 'statd.for';
C-        _include grid;
C: *
C---- '-->Entering Node %_newphi_vars:'
      character*(*) rotlaw
      Integer im,ip,k,KMax,kk,KA,KB,K0,L,IL,IRN,ICN,ID,IDD,IDISP,IFAIL,
     *	IKEEP,I00,ITQ,IVECT,IW,JVECT,LICN,LIRN,MTYPE,NIN,NOUT,nnzero,i
     *    rcylb,ir0,ic0,kbphIm,kbphIp,n,NZ
      doubleprecision A,b,b0,bim,bip,bjm,bjp,bK0,bKA,bKB,bphim,bphip,
     *	ff,a_1,b_1,C,c_1,a_2,b_2,c_2,dC_dphi0,dC_dphiA,dC_dpsiA,
     *    drho_dC,drho_dH,drho_dphi,drho_dphiu,drho_dpsi,EOS,phircyl,
     *    dphiu_dC,dphiu_dphi0,psi,dpsi0_dphi0,dpsi0_dphiA,dpsi0_dphiB,
     *    dpsi0_dphim,dpsi0_dphip,dpsi0A_dphi0,dpsi0A_dphiA,dpsi0A_dphiB
     *    ,dpsi_dphiA,dpsi_dphiB,dpsi_dpsi0,dpsiA_dphiA,dpsiA_dphiB,
     *    dpsiA_dpsi0,dpsiA_dpsi0A
      doubleprecision Asave,Bff,chi,chiNorm,chi_jmscf,chi_solid,chimax,
     *	chiwork,ERREST,ERROLD,ETA,FSAVE,Fpsi,Fpsi0,powrsmart,potential
     *    ,potential_min,RESID,RPMIN,rcylm,rcylp,theta1,theta2,thetaa,
     *    THR,Upivot,W,Xsolution,XMOD,xm,xp
      Parameter(KMax=(IMax+1)*(JMax+1))
      Parameter(nnzero=20*KMax)
      Parameter(LICN=10*nnzero,LIRN=10*nnzero/2)
C-        Parameter(LICN=10*nnzero,LIRN=10*nnzero);
      Parameter(NIN=5,NOUT=1)
      Dimension
C- r(0:IMax),alpha(0:IMax),theta(0:JMax),
     *   a_1(0:IMax),b_1(0:IMax),c_1(0:IMax),a_2(0:JMax),b_2(0:JMax),
     *   c_2(0:JMax),b(KMax),bim(KMax),bip(KMax),bjm(KMax),bjp(KMax),
     *   bK0(KMax),bKA(KMax),bKB(KMax),b0(0:JMax),bphim(KMax),
     *   bphip(KMax),ff(KMax)
      Dimension kbphIm(KMax),kbphIp(KMax)
      Dimension ir0(nnzero),ic0(nnzero)
      Dimension A(LICN),W(KMax),ICN(LICN),IDISP(11),IKEEP(KMax,5),
     *	IRN(LIRN),IW(KMax,8),IVECT(nnzero),JVECT(nnzero)
      Dimension Asave(nnzero),Fsave(Kmax),Bff(KMax),Xsolution(Kmax)
      Integer MAXIT
      Doubleprecision ANOISE
      Parameter(ANOISE=1.D-11)
      Parameter(MAXIT=30)
      Logical GROW,LBLOCK,ABORT(4)
      Logical LogInit,LogRapid
C-        Common /phi/phi(0:IMax,0:JMax);
C-        Common /delta/delta;
C-        Common /pi/pi;
C---- '<--Leaving  Node %_newphi_vars:'
C       dr=one/IA;
C       dalpha=one/r(IA-1)-one/r(IA);*
      KA=kk(IA,JMax)
      KB=kk(IB,JB)
      K0=kk(I0,JMax)
      do i=1,KMax
      b(i)=zero
      bim(i)=zero
      bip(i)=zero
      bjm(i)=zero
      bjp(i)=zero
      ff(i)=zero
      enddo
      do i=1,IA-1
      a_1(i)=r(i)**(one-three*delta)/(delta*dr)**2*
     *	(half*(r(i-1)+r(i)))**(one+delta)
      c_1(i)=r(i)**(one-three*delta)/(delta*dr)**2*
     *	(half*(r(i)+r(i+1)))**(one+delta)
      b_1(i)=-a_1(i)-c_1(i)
      enddo
      do l=1,IMax-IA
      i=IMax-l
      a_1(i)=alpha(l)**(one+three*delta)/(delta*dalpha)**2*
     *	(half*(alpha(l)+alpha(l+1)))**(one-delta)
      c_1(i)=alpha(l)**(one+three*delta)/(delta*dalpha)**2*
     *	(half*(alpha(l-1)+alpha(l)))**(one-delta)
      b_1(i)=-a_1(i)-c_1(i)
      enddo
      do j=0,JMax
      if((j.eq.0).or.(j.eq.JMax))then
      if(j.eq.0)then
      c_2(j)=two**2/dtheta**2
      b_2(j)=-c_2(j)
      else
      a_2(j)=two/dtheta**2
      b_2(j)=-a_2(j)
      endif
      else
      a_2(j)=one/(sin(theta(j))*dtheta**2)*
     *	sin(half*(theta(j-1)+theta(j)))
      c_2(j)=one/(sin(theta(j))*dtheta**2)*
     *	sin(half*(theta(j)+theta(j+1)))
      b_2(j)=-a_2(j)-c_2(j)
      endif
      enddo
C-     === approximation of Laplace-operator ==>
C-     === internal space ==>
      do i=1,IMax-1
      do j=0,JMax
      	k=kk(i,j)
      	b(k)=b_1(i)+b_2(j)/r(i)**(2*delta)
      	bim(k)=a_1(i)
      bip(k)=c_1(i)
      	if(j.ne.0)then
      bjm(k)=a_2(j)/r(i)**(2*delta)
      endif
      	if(j.ne.JMax)then
      		bjp(k)=c_2(j)/r(i)**(2*delta)
      endif
      enddo
      enddo
C-     <== internal space ===
C-     === r=0 ==>
      b(kk(0,JMax))=-6d0/dr**(2*delta)
      do j=0,JMax
      theta1=half*(theta(max(0,j-1))+theta(j))
      theta2=half*(theta(j)+theta(min(JMax,j+1)))
      	b0(j)=6d0/dr**(2*delta)*(cos(theta1)-cos(theta2))
      enddo
      do j=0,JMax-1
      	k=kk(0,j)
      	b(k)=one
      	bjp(k)=-one
      enddo
C-     <== r=0 ===
C-     === alpha=0 ==>
      do j=0,JMax
      	b(kk(IMax,j))=one
      enddo
C-     <== alpha=0 ===
C-     <== approximation of Laplace-operator ===
      do i=0,IA
      do j=0,JMax
      	k=kk(i,j)
      	if(k.ge.kk(0,JMax))then
      	psi=Fpsi(powrsmart(r(i),delta)*sin(theta(j)),rotlaw)
      rho=EOS(phiu*(C-phi(i,j)-psi),1)/rho0
C- key=1 see Eq. 18
C-               write(*,'(a,2i5,1p,4g12.3)')' i,j, C, phi,psi,rho=',i,j,C,phi(i,j), psi,rho;
C-             pause;
      drho_dH=EOS(phiu*(C-phi(i,j)-psi),2)
C- key=2
      drho_dphiu=drho_dH/rho0*(C-phi(i,j)-psi)
C- see Eq.9
      drho_dC=drho_dH/rho0*phiu
      drho_dphi=-drho_dH/rho0*phiu
      drho_dpsi=-drho_dH/rho0*phiu
      dphiu_dC=-EOS(rho0,3)/(C-phi(I0,JMax))**2
      dphiu_dphi0=EOS(rho0,3)/(C-phi(I0,JMax))**2
      dC_dphi0=(-EOS(zero,3)/EOS(rho0,3))/(one-EOS(zero,3)/EOS(rho0,3))
C- see Eq. 17
      dC_dphiA=one/(one-EOS(zero,3)/EOS(rho0,3))
      dC_dpsiA=one/(one-EOS(zero,3)/EOS(rho0,3))
      chiNorm=1.d0/(chi(powrsmart(r(IA),delta)*sin(theta(JMax)),
     *	rotlaw(1:len_trim(rotlaw)))-chi(powrsmart(r(IB),delta)*
     *    sin(theta(JB)),rotlaw(1:len_trim(rotlaw))))
      dpsi_dphiA=-chi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))*chiNorm
C- as eq.14, 16
      dpsi_dphiB=chi(powrsmart(r(i),delta)*sin(theta(j)),
     *	rotlaw(1:len_trim(rotlaw)))*chiNorm
      dpsi_dpsi0=phi(IB,JB)-phi(IA,JMax)
      dpsiA_dphiA=-chi(one,rotlaw)*chiNorm
      dpsiA_dphiB=chi(one,rotlaw)*chiNorm
      dpsiA_dpsi0=phi(IB,JB)-phi(IA,JMax)
      dpsiA_dpsi0A=dpsiA_dpsi0
      if(LogRapid)then
      if(JB.eq.0)then
      ircylb=0
      else
      ircylb=IB
      endif
      rcyl=powrsmart(r(i),delta)*sin(theta(j))
      im=powrsmart(rcyl,one/delta)/dr
      ip=im+1
      rcylm=im*dr
      rcylp=ip*dr
      xm=(rcylp-powrsmart(rcyl,one/delta))/dr
      xp=(powrsmart(rcyl,one/delta)-rcylm)/dr
      phircyl=xm*phi(im,JMax)+xp*phi(ip,JMax)
      dpsi0_dphi0=-one/(phi(IA,JMax)-phi(ircylb,JMax))
      dpsi0A_dphi0=-one/(phi(IA,JMax)-phi(ircylb,JMax))
      dpsi0_dphiA=-(phircyl-phi(I0,JMax))/(phi(IA,JMax)-
     *	phi(ircylb,JMax))**2
      dpsi0A_dphiA=-(phi(IA,JMax)-phi(I0,JMax))/(phi(IA,JMax)-
     *	phi(ircylb,JMax))**2+one/(phi(IA,JMax)-phi(ircylb,JMax))
      if(JB.eq.0)then
      	dpsi0_dphiB=zero
      dpsi0_dphi0=dpsi0_dphi0+(phircyl-phi(I0,JMax))/(phi(IA,JMax)-
     *	phi(ircylb,JMax))**2
      	dpsi0A_dphiB=zero
      dpsi0A_dphi0=dpsi0A_dphi0+(phi(IA,JMax)-phi(I0,JMax))/
     *	(phi(IA,JMax)-phi(ircylb,JMax))**2
      else
      	dpsi0_dphiB=(phircyl-phi(I0,JMax))/(phi(IA,JMax)-
     *	phi(ircylb,JMax))**2
      	dpsi0A_dphiB=(phi(IA,JMax)-phi(I0,JMax))/(phi(IA,JMax)-
     *	phi(ircylb,JMax))**2
      endif
      	dpsi0_dphim=xm/(phi(IA,JMax)-phi(ircylb,JMax))
      	dpsi0_dphip=xp/(phi(IA,JMax)-phi(ircylb,JMax))
      else
      dpsi0_dphi0=zero
      dpsi0A_dphi0=zero
      dpsi0_dphiA=zero
      dpsi0A_dphiA=zero
      dpsi0_dphiB=zero
      dpsi0A_dphiB=zero
      	dpsi0_dphim=zero
      	dpsi0_dphip=zero
      endif
      ff(k)=4d0*pi*(rho-drho_dphiu*(dphiu_dphi0*phi(I0,JMax)+dphiu_dC*
     *	(dC_dphi0*phi(I0,JMax)+dC_dphiA*phi(IA,JMax)+dC_dpsiA*(
     *    dpsiA_dphiA*phi(IA,JMax)+dpsiA_dphiB*phi(IB,JB)+dpsiA_dpsi0A*
     *   (dpsi0A_dphi0*phi(I0,JMax)+dpsi0A_dphiA*phi(IA,JMax)+
     *    dpsi0A_dphiB*phi(IB,JB)))))-drho_dC*(dC_dphi0*phi(I0,JMax)+
     *    dC_dphiA*phi(IA,JMax)+dC_dpsiA*(dpsiA_dphiA*phi(IA,JMax)+
     *    dpsiA_dphiB*phi(IB,JB)+dpsiA_dpsi0A*(dpsi0A_dphi0*phi(I0,JMax)
     *   +dpsi0A_dphiA*phi(IA,JMax)+dpsi0A_dphiB*phi(IB,JB))))-drho_dphi
     *   *phi(i,j)-drho_dpsi*(dpsi_dphiA*phi(IA,JMax)+dpsi_dphiB*
     *   phi(IB,JB)+dpsi_dpsi0*(dpsi0_dphi0*phi(I0,JMax)+dpsi0_dphiA*
     *   phi(IA,JMax)+dpsi0_dphiB*phi(IB,JB)+dpsi0_dphim*phi(im,JMax)+
     *   dpsi0_dphip*phi(ip,JMax))))
      b(k)=b(k)+4d0*pi*(-drho_dphi)
      bK0(k)=+4d0*pi*(-drho_dphiu*(dphiu_dphi0+dphiu_dC*(dC_dphi0+
     *	dC_dpsiA*(+dpsiA_dpsi0A*(dpsi0A_dphi0))))-drho_dC*(dC_dphi0+
     *    dC_dpsiA*(+dpsiA_dpsi0A*(dpsi0A_dphi0)))-drho_dpsi*(
     *    +dpsi_dpsi0*(dpsi0_dphi0)))
      bKA(k)=+4d0*pi*(-drho_dphiu*(+dphiu_dC*(+dC_dphiA+dC_dpsiA*
     *	(dpsiA_dphiA+dpsiA_dpsi0A*(+dpsi0A_dphiA))))-drho_dC*(
     *    +dC_dphiA+dC_dpsiA*(dpsiA_dphiA+dpsiA_dpsi0A*(+dpsi0A_dphiA)))
     *    -drho_dpsi*(dpsi_dphiA+dpsi_dpsi0*(+dpsi0_dphiA)))
      bKB(k)=+4d0*pi*(-drho_dphiu*(+dphiu_dC*(+dC_dpsiA*(+dpsiA_dphiB+
     *	dpsiA_dpsi0A*(+dpsi0A_dphiB))))-drho_dC*(+dC_dpsiA*(
     *   +dpsiA_dphiB+dpsiA_dpsi0A*(+dpsi0A_dphiB)))-drho_dpsi*(
     *   +dpsi_dphiB+dpsi_dpsi0*(+dpsi0_dphiB)))
      bphim(k)=+4d0*pi*(-drho_dpsi*(+dpsi_dpsi0*(+dpsi0_dphim)))
      kbphim(k)=kk(im,JMax)
      bphip(k)=+4d0*pi*(-drho_dpsi*(+dpsi_dpsi0*(+dpsi0_dphip)))
      	kbphip(k)=kk(ip,JMax)
      	endif
      enddo
      enddo
C- <== internal space ===
C-     === forming initial data for sparse matrix solver ==>
      NZ=0
      do k=1,KMax
      i=(k-1)/(JMax+1)
      	j=k-1-i*(JMax+1)
      NZ=NZ+1
      	A(NZ)=b(k)
C- main diagonal
      	ir0(NZ)=k
      	ic0(NZ)=k
      if(i.ne.0)then
      NZ=NZ+1
      	A(NZ)=bim(k)
C- parall. to main diagonal
      	ir0(NZ)=k
      	ic0(NZ)=kk(i-1,j)
      endif
      if(i.ne.IMax)then
      NZ=NZ+1
      	A(NZ)=bip(k)
C- parall. to main diagonal
      	ir0(NZ)=k
      	ic0(NZ)=kk(i+1,j)
      endif
      if(j.ne.0)then
      NZ=NZ+1
      	A(NZ)=bjm(k)
C- parall. to main diagonal
      	ir0(NZ)=k
      	ic0(NZ)=kk(i,j-1)
      endif
      if(j.ne.JMax)then
      NZ=NZ+1
      	A(NZ)=bjp(k)
C- parall. to main diagonal
      	ir0(NZ)=k
      	ic0(NZ)=kk(i,j+1)
      endif
      	if((i.le.IA).and.(k.ge.kk(0,JMax)))then
      NZ=NZ+1
      	A(NZ)=bK0(k)
C- point 0
      	ir0(NZ)=k
      	ic0(NZ)=K0
      NZ=NZ+1
      	A(NZ)=bKA(k)
C- point A
      	ir0(NZ)=k
      	ic0(NZ)=KA
      NZ=NZ+1
      	A(NZ)=bKB(k)
C- point B
      	ir0(NZ)=k
      	ic0(NZ)=KB
      if(LogRapid)then
C- alpha-law
      NZ=NZ+1
      A(NZ)=bphim(k)
      ir0(NZ)=k
      ic0(NZ)=kbphim(k)
      NZ=NZ+1
      A(NZ)=bphip(k)
      ir0(NZ)=k
      ic0(NZ)=kbphip(k)
      endif
      	endif
      enddo
      do j=0,JMax
      NZ=NZ+1
      A(NZ)=b0(j)
C- \nabla^2\Phi at r->0
      	ir0(NZ)=kk(0,JMax)
      	ic0(NZ)=kk(1,j)
      enddo
      do n=1,NZ
      if(LogInit)then
      	IRN(n)=ir0(n)
      	ICN(n)=ic0(n)
      else
      	IVECT(n)=ir0(n)
      JVECT(n)=ic0(n)
      endif
      Asave(n)=A(n)
      enddo
      Upivot=0.1d0
C-        Upivot=0.01d0;
      LBLOCK=.TRUE.
      GROW=.TRUE.
      ABORT(1)=.TRUE.
      ABORT(2)=.TRUE.
      ABORT(3)=.FALSE.
      ABORT(4)=.FALSE.
C- IFAIL=110
      IFAIL=010
C-     <== end forming initial data for sparse matrix solver ==>
      CALL X04ABF(1,NOUT)
      if(LogInit)then
C-**          CALL F01BRF(KMax,NZ,A,LICN,IRN,LIRN,ICN,Upivot,IKEEP,IW,W,LBLOCK,
C-**                     GROW,ABORT,IDISP,IFAIL); -- old NAG sparse solver
      THR=1d-14
CCC      THR=1d-24
      CALL M28Y12(KMax,NZ,A,LICN,IRN,LIRN,ICN,Upivot,IKEEP,IW,W,LBLOCK,
     *	GROW,ABORT,IDISP,IFAIL,THR)
      WRITE(*,99996)IDISP(2),IDISP(6),IDISP(7),IDISP(3),IDISP(4)
      IF(LBLOCK)WRITE(*,99994)(IDISP(I),I=8,10)
99996 FORMAT(' NUMBER OF NON-ZEROS IN DECOMPOSITION',9X,'=',I8/'
     *    MINIMUM SIZE OF ARRAY IRN',20X,'=',I8/' MINIMUM SIZE OF ARRAYS
     *    A AND ICN',13X,'=',I8/' NUMBER OF COMPRESSES ON IRN (
     *    IDISP(3))',7X,'=',I8/' NUMBER OF COMPRESSES ON A AND ICN
     *    (IDISP(4)) =',I8)
99994 FORMAT(' STRUCTURAL RANK',16X,'=',I7/' NUMBER OF DIAGONAL BLOCKS'
     *   ,6X,'=',I7/' SIZE OF LARGEST DIAGONAL BLOCK =',I7)
      LogInit=.false.
C-          LogInit=.true. ! only M28BYS (25 October)
C-          NM28BYS=0
      else
      ETA=1d-4
C- IFAIL=110
      IFAIL=010
C-**          CALL F01BSF(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,
C-**                     ETA,RPMIN,ABORT(4),IDISP,IFAIL); -- old NAG sparse solver
      CALL M28BYS(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,ETA,
     *	RPMIN,ABORT(4),IDISP,IFAIL)
C-          NM28BYS=NM28BYS+1
C-          if(NM28BYS.ge.2) LogInit=.true.
      endif
      MTYPE=1
C:  iterative refinement by M28CYS or M28CYN *
C---- '-->Entering Node %_newphi_sol:'
C- iterative refinement by M28CYS/N
      DO 09999 IDD=1,Kmax
      Bff(IDD)=ff(IDD)
C- R-H SIDE IN Bff, ff saves RHS
      Xsolution(IDD)=0.D0
C- clear Xsolution
09999 CONTINUE
      ITQ=0
C-     write(*,'(A)' )' The run of estimated error ERREST in iterations:';
09996 CONTINUE
      CALL M28CYN(Kmax,A,LICN,ICN,IKEEP,Bff,W,IW,1,IDISP,RESID,Ifail)
C-CALL M28CYS(KMax,A,LICN,ICN,IKEEP,ff,W,MTYPE,IDISP,RESID);
C-CALL M28CYS(Kmax,A,LICN,ICN,IKEEP,Bff,W,MTYPE,IDISP,RESID);
C- solution is in Bff
      If(Ifail.NE.0)then
      write(*,'(A,I2)')' Ifail in M28CYN =',Ifail
      endif
      ERROLD=ERREST
      ERREST=0.D0
      XMOD=0.D0
      DO 09993 ID=1,Kmax
      Xsolution(ID)=Xsolution(ID)+Bff(ID)
      XMOD=XMOD+ABS(Xsolution(ID))
C- NORMA
      ERREST=ERREST+ABS(Bff(ID))
      FSAVE(ID)=0.D0
09993 CONTINUE
      DO 09990 K=1,NZ
      FSAVE(IR0(K))=FSAVE(IR0(K))+Xsolution(IC0(K))*ASAVE(K)
09990 CONTINUE
      DO 09987 IL=1,KMax
      Bff(IL)=ff(IL)-FSAVE(IL)
C- Polish strategy
09987 CONTINUE
C-       write(*,'(A,I2,A,1P,E10.3)' )' It=',ITQ,'    ERREST=',ERREST;
      ITQ=ITQ+1
      IF(.NOT.(ERREST.LE.ANOISE*XMOD.or.(ITQ.GT.2.AND.ERREST.GT.ERROLD)
     *	.or.ITQ.GT.MAXIT))GOTO 09996
      DO 09984 IL=1,Kmax
      ff(IL)=Xsolution(IL)
C- Polish strategy
09984 CONTINUE
C---- '<--Leaving  Node %_newphi_sol:'
      do k=1,KMax
      i=(k-1)/(JMax+1)
      	j=k-1-i*(JMax+1)
      phi(i,j)=ff(k)
      enddo
      if(JB.eq.JMax)then
C- toroidal configurations
C- let new  I0
      I00=I0
      potential_min=phi(I0,JMax)+Fpsi(powrsmart(r(I0),delta),rotlaw)
      do i=IB+1,IA-1
      potential=phi(i,JMax)+Fpsi(powrsmart(r(i),delta),rotlaw)
      if(potential.LT.potential_min)then
      I0=i
      potential_min=potential
      endif
      enddo
      if(I00.ne.I0)LogInit=.false.
      endif
      C=(phi(IA,JMax)+Fpsi(one,rotlaw(1:len_trim(rotlaw)))-phi(I0,JMax)*
     *	EOS(zero,3)/EOS(rho0,3)
C- EOS(*,3) means enthalpy H
     *)/(one-EOS(zero,3)/EOS(rho0,3))
C- eq.17 in AA290(1994)674
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
C- eq.18 assuming A at r==0, what if r>0?
      return
      end
C---- '<--Leaving  Node %_newphi:'
C:  integer Function kk(i,j); -- simply finds
C       kk=1+i*(JMax+1)+j;     *
C---- '-->Entering Node %_kk:'
      integer Function kk(i,j)
	use cylinderGrid, only: IA,JMax,IMax
      Implicit none
C-*****************************************************************
C-       include 'statd.for';
      integer,Parameter::irho0max=5
C- 50 standard  -- 5 for test runs
CCC      integer::IA,JMax,IMax
C-      Parameter(IA  =600,JMax=10);
C-      Parameter(IA  =600,JMax=5); -- for non-rotating
C-      Parameter(IA  =400,JMax=10); -- for non-rotating
C-      Parameter(IA  =600,JMax=10);
C-      Parameter(IA  =200,JMax=100); -- for rotating
C-      Parameter(IA  =200,JMax=5); -- for non-rotating
C      Parameter(IA=200,JMax=100);
C-       Parameter(IA  =200,JMax=50);
C-       Parameter(IA  =40,JMax=5);
C      Parameter(IA=5000,JMax=80);
C      Parameter(IA=150,JMax=80)
CCC      Parameter(IA=200,JMax=100)
C- for fast crude tests
CCC      Parameter(IMax=2*IA-1)
      integer::i,j
      kk=1+i*(JMax+1)+j
      return
      end
C---- '<--Leaving  Node %_kk:'
C:        Fpsi=(phi(IB,JB)-phi(IA,JMax))*chi(rcyl); -- almost eq.14, but C_\Psi
C             -- is replaced by phi(IB,JB)-phi(IA,JMax) -- можно ли? наверное только для alpha-law
C        *
C---- '-->Entering Node %_fpsi:'
      doubleprecision Function Fpsi(rcylw,rotlaw)
C- centrifugal potential
C-*****************************************************************
      use accur
      use cylinderGrid
      use utilities_module
      Implicit none
C-        _include accur;
C-        include 'accur.for';
C-       include 'statd.for';
C-        _include grid;
C-        Integer IB,JB;
      doubleprecision chi,rcylw,powrsmart
      character*(*)rotlaw
C-        Common /phi/phi(0:IMax,0:JMax);
C-        Common /IB/IB/JB/JB;
C-       write(*,'(a,4i5,1pg12.3,2x,a)')'Fpsi: IB,JB,IA,JMax,rcylw,rotlaw =',IB,JB,IA,JMax,rcylw,rotlaw;
      if(rotlaw(1:len_trim(rotlaw)).EQ.'alpha')then
      Fpsi=(phi(IB,JB)-phi(IA,JMax))*chi(rcylw,
     *	rotlaw(1:len_trim(rotlaw)))
C- almost eq.14,16  but C_\Psi
C- is replaced by phi(IB,JB)-phi(IA,JMax) as Eq.55 -- for alpha-law
      else
      Fpsi=(phi(IB,JB)-phi(IA,JMax))*chi(rcylw,rotlaw(1:len_trim(rotlaw
     *	)))/(chi(powrsmart(r(IA),delta)*sin(theta(JMax)),rotlaw(1:len
     *_trim(rotlaw)))-chi(powrsmart(r(IB),delta)*sin(theta(JB)),
     * rotlaw(1:len_trim(rotlaw))))
C- as eq.14, 16
      endif
      Cpsi =(phi(IB,JB)-phi(IA,JMax))/
     * (chi(powrsmart(r(IA),delta)*sin(theta(JMax)),
     * rotlaw(1:len_trim(rotlaw)))-chi(powrsmart(r(IB),delta)*
     * sin(theta(JB)),rotlaw(1:len_trim(rotlaw))))

C-       write(*,'(a,1p4g12.3)')'Fpsi phi_B phi_A,chi =',Fpsi,phi(IB,JB),phi(IA,JMax),chi(rcyl,rotlaw);
C-       call mypause;
      return
      end
C---- '<--Leaving  Node %_fpsi:'
C: rotation potential:
C e.g. eq. 15
C $\hat\Psi^0(\hat\varpi)$
Cis a prescribed function, for example (Eriguchi \& M\"uller 1985a),
C\begin{equation}
C    \hat\Psi^0(\hat\varpi)
C        =\cases{-\hat\varpi^2/2 -- here better plus
C                & rigid rotation,
C                \cr
C                -\frac{1}{2} -- here better plus
C                 \ln(d^2+\hat\varpi^2)
C                & $v$-constant rotation,
C                \cr
C                -1/[2(d^2+\hat\varpi^2)]
C                & $j$-constant rotation;  -- sign in chi may change
C         } \label{eq:rolaw}
C\end{equation}
Chere $d$ is some constant.
Cor eq.50
C   *
C---- '-->Entering Node %_chi:'
      doubleprecision Function chi(rcylw,rotlaw)
C- \hat\Psi^0
C-*****************************************************************
C- chi(r(I0))=zero
C- chi(r(IA))-chi(powrsmart(r(IB),delta)*sin(theta(JB)))=one -- this normalization
C- is not necessary when Fpsi divides by chi(r(IA))-chi(powrsmart(r(IB),delta)*sin(theta(JB)))
      use accur
      use cylinderGrid
      Implicit none
C-        _include accur;
C-        include 'accur.for';
C-       include 'statd.for';
C-        _include grid;
      Integer ircylb,im,ip
      doubleprecision cjw,chimax,chiwork,powrsmart,rcylw,rcylm,rcylp
      doubleprecision,parameter::dsquared=0.01d0
C- from Hachisu 1986
      doubleprecision myrotlaw,tau,beta
	doubleprecision omegalaw,savedlaw,jfixedlaw
	data tau/0.7D0/beta/0.9D0/
      character*(*)rotlaw

C-        Common /phi/phi(0:IMax,0:JMax);
C-        Common /delta/delta;
C-        Common /IB/IB/JB/JB/I0/I0;
C-        Logical firstJmscf/.true/;
C-        save firstJmscf;
      if(JB.eq.0)then
      ircylb=0
      else
      ircylb=IB
      endif
      if(rotlaw(1:len_trim(rotlaw)).EQ.'rigid')then
      if(rcylw.le.one)then
      chi=+half*rcylw**2
C- changing sign, because no C_\psi in Eq.16 is defined as in mrot*
C- test another constant
C-                  chi=+1d3*rcylw**2; -- changing sign, because no C_\psi in Eq.16 is defined as in mrot*
      else
      chi=+half
C-                  chi=+1d3;
      endif
C-        write(*,*)' in chi, chi, rotlaw:', chi,rotlaw,'  ',rotlaw(1:len_trim(rotlaw));
      elseif(rotlaw(1:len_trim(rotlaw)).EQ.'jmscf')then
C- not debugged yet
      if(LinitMcyl)then
C- repeating rigid
      if(rcylw.le.one)then
      chi=+half*rcylw**2
      else
      chi=+half
      endif
      else
C: integrate j(mcyl)=2.5-2.5*(1-mcyl)**(2/3) to get chi and Fpsi  *
C---- '-->Entering Node %_chi_jmscf:'
C- j(mcyl)=2.5-2.5*(1-mcyl)**(2/3)
C- chi must be \int_0^rcyl ds j**2/s**3 -- use GillMiller for that?
C- here j=j(mcyl) with mcyl(rcyl)= \int_0^rcyl ds sigma(s) s / total mass
C-  CMcyl(0:NrDAT),CJcyl(0:NrDAT),CMcylg(0:NrDAT),CJcylg(0:NrDAT), -- suffix g for grid values
C-
      do icyl=1,NrDat
C- always begin from the outer radius to normalize
C-              cjw=2.5d0-2.5d0*(1.d0-mcyl(icyl))**(2.d0/3.d0);
      cjw=1.d0-(1.d0-mcyl(icyl))**(2.d0/3.d0)
      chimax=chimax+cjw**2/Rcyld(icyl)**3
      enddo
      open(101,file='testJm.res')
      write(101,'(a)')'#     Rcyl           chi_solid      chi_jmscf'
      chiwork=0.d0
      do icyl=1,NrDat
      cjw=1.d0-(1.d0-mcyl(icyl))**(2.d0/3.d0)
      chiwork=chiwork+cjw**2/Rcyld(icyl)**3
      chi=half*chiwork/chimax
      write(101,'(1p,3g15.3)')Rcyld(icyl),half*Rcyld(icyl)**2,chi
      enddo
      stop 'testJm.res written'
      if(rcylw.LT.Rcyld(1))then
      chi=0.d0
      else
      chiwork=0.d0
      icyl=1
      do while(Rcyld(min(icyl,NrDAT)).LE.rcylw.and.icyl.LE.NrDAT)
C-            do while  (Rcyld(icyl) <= rcylw .and. icyl<=NrDAT);
C-              cjw=2.5d0-2.5d0*(1.d0-mcyl(icyl))**(2.d0/3.d0);
      cjw=1.d0-(1.d0-mcyl(icyl))**(2.d0/3.d0)
      chiwork=chiwork+cjw**2/Rcyld(icyl)**3
C             if(icyl==1 .or. icyl==NrDAT)then;
C               write(*,'(a,i5,1p,5g15.5)')' icyl,mcyl(icyl),cjw,chi,rcylw,Rcyld(icyl)=',icyl,mcyl(icyl),cjw,chiwork,rcylw,Rcyld(icyl);
C             endif;*
      icyl=icyl+1
      enddo
C- exit when Rcyld(i) >= rcylw
      endif
      if(rcylw.le.one)then
      chi=half*chiwork/chimax
      else
      chi=+half
      endif
C-           write(*,'(a,1p,5g15.5)')'half*rcylw**2,chi', half*rcylw**2,chi;
C-           pause;
C
C          write(*,'(a,i5,1p,5g15.5)')' icyl,mcyl(icyl),cjw,chi,rcylw,Rcyld(icyl)=',
C                 icyl,mcyl(icyl),cjw,chi,rcylw,Rcyld(icyl);
C          pause;
C          *
C---- '<--Leaving  Node %_chi_jmscf:'
      endif
      elseif(rotlaw(1:len_trim(rotlaw)).EQ.'vconst')then
C- from Hachisu 1986 for d^2=0.01
      if(rcylw.le.one)then
      chi=+half*log(dsquared+rcylw**2)
C- changing sign, because no C_\psi in Eq.16 is defined as in mrot*
      else
      chi=+half*log(dsquared+one)
      endif
      elseif(rotlaw(1:len_trim(rotlaw)).EQ.'jconst')then
C- from Hachisu 1986 for d^2=0.01
      if(rcylw.le.one)then
      chi=-half*one/(dsquared+rcylw**2)
C- here minus sign, if we want growing chi
      else
      chi=-half*one/(dsquared+one)
      endif
	elseif(rotlaw(1:len_trim(rotlaw)).eq.'collap')then
	chi=myrotlaw(rcylw,tau,beta)
	elseif(rotlaw(1:len_trim(rotlaw)).eq.'omega')then
	chi=omegalaw(rcylw)
	elseif(rotlaw(1:len_trim(rotlaw)).eq.'saved')then
	chi=savedlaw(rcylw)
      elseif(rotlaw(1:len_trim(rotlaw)).eq.'jfixed')then
	chi=jfixedlaw(rcylw)
      else
C- alpha-law Eq.50
      if(rcylw.lt.one)then
      dr=one/dble(IA)
      im=powrsmart(rcylw,one/delta)/dr
      ip=im+1
      rcylm=im*dr
      rcylp=ip*dr
C-                   Fpsi0:
      chi=((rcylp-powrsmart(rcylw,one/delta))/dr*phi(im,JMax)+(powrsmart
     *	(rcylw,one/delta)-rcylm)/dr*phi(ip,JMax)-phi(I0,JMax))/(phi
     *(IA,JMax)-phi(ircylb,JMax))
      else
C-                  Fpsi0:
C-                  write(*,*)' IA,JMax,I0,ircylb=', IA,JMax,I0,ircylb;
      chi=(phi(IA,JMax)-phi(I0,JMax))/(phi(IA,JMax)-phi(ircylb,JMax))
      endif
      endif
C-        write(*,*)' chi,rcylw=',chi,rcylw;
C-        pause;
C-        stop;
      return
      end
C---- '<--Leaving  Node %_chi:'
C: numerical derivative of Psi over cyl.radius
C        double precision Function Fdpsi_drcyl(rcyl);
C        drcyl=1d-5;
C           rcylm=max(zero,rcyl-half*drcyl);
C           rcylp=min(one,rcyl+half*drcyl);
C           Fdpsi_drcyl=(Fpsi(rcylp,rotlaw)-Fpsi(rcylm,rotlaw))/(rcylp-rcylm);
C           -- better to do analytcs for simple rotation laws
C  *
C---- '-->Entering Node %_fdpsi:'
C:   *
C---- '-->Entering Node %_fdpsi_drcyl:'
      doubleprecision Function Fdpsi_drcyl(rcylw,rotlaw)
C-*****************************************************************
      use accur
      use cylinderGrid
      Implicit none
C-        _include accur;
C-        include 'accur.for';
C-       include 'statd.for';
C-        _include grid;
      doubleprecision drw,drcylw,Fpsi,rcylw,rcylm,rcylp
      character*(*)rotlaw
      drw=one/IA
C- put to common or module, or arg
      if(rcylw.le.one)then ! Mine change 18.07.2018 .lt. --> .le.
      drcylw=1.d-5*drw
      rcylm=max(zero,rcylw-half*drcylw)
      rcylp=min(one,rcylw+half*drcylw)
      Fdpsi_drcyl=(Fpsi(rcylp,rotlaw)-Fpsi(rcylm,rotlaw))/(rcylp-rcylm)
      else
      Fdpsi_drcyl=zero
      endif
      return
      end
C---- '<--Leaving  Node %_fdpsi_drcyl:'
C---- '<--Leaving  Node %_fdpsi:'
C: *
C---- '-->Entering Node %_powrsmart:'
C-       double precision Function degree(x,alpha); -- Aksenov's name for this function
      doubleprecision Function powrsmart(x,alpha)
C-*****************************************************************
      use accur
      doubleprecision alpha,x
C-        include 'accur.for';
      if(alpha.eq.zero)then
      powrsmart=one
      else
      if(x.eq.zero)then
      powrsmart=zero
      else
      powrsmart=x**alpha
      endif
      endif
      return
      end
C---- '<--Leaving  Node %_powrsmart:'
C:   *
C---- '-->Entering Node %_integr:'
      Subroutine Integrals(CM,V,CJ,CI,T,W,Eth,Pip,C,omegaA,rotlaw)
C-*****************************************************************
      use accur
      use cylinderGrid
      Implicit none
C-       _include accur;
C-        include 'accur.for';
C-       include 'statd.for';
C-        _include grid;
C-        integer i,j,IB,JB;
      doubleprecision CM,V,dV,CJ,CI,T,W,Eth,P,Pip,ra,thetaa,C,psi,Fpsi,
     *	Fdpsi_drcyl,omegaA
      doubleprecision c1,powrsmart,EOS
C-        Dimension r(0:IMax),theta(0:JMax),
C-                  phi(0:IMax,0:JMax);
C-        Common /pi/pi/IB/IB/JB/JB;
C-        Common /IB/IB/JB/JB;
C-        Common /delta/delta;
C-        Common /gamma/gamma,CN,CK;
      character*(*)rotlaw
      CM=zero
C- mass
      V=zero
C- volume
      CJ=zero
C- angular momentum
      CI=zero
C- moment of inertia
      T=zero
C- kinetic energy
      W=zero
C- grav.pot. energy
      Pip=zero
C- pressure integral for virial theorem
      Eth=zero
C- thermal energy
      do i=1,IA
      do j=1,JMax
      ra=half*(r(i-1)+r(i))
C- different from eq.43 in AA290(1994)674 (no max/min)
C- cf. mrot.trf %_integrals where it is closer to eq.43
      thetaa=half*(theta(j-1)+theta(j))
      phia=(phi(i-1,j-1)+phi(i-1,j)+phi(i,j-1)+phi(i,j))/four
      psi=Fpsi(powrsmart(ra,delta)*sin(thetaa),
     *	rotlaw(1:len_trim(rotlaw)))
      rho=EOS(phiu*(C-phia-psi),1)/rho0
      c1=-Fdpsi_drcyl(powrsmart(ra,delta)*sin(thetaa),
     *	rotlaw(1:len_trim(rotlaw)))*powrsmart(ra,delta)*sin(thetaa)
C- it is minus, now the sign is opposite in chi -- check
      	vphi=sqrt(max(zero,c1))
C-               write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,c1,vphi;
C-               pause;
      dV=(powrsmart(r(i),delta)**3-powrsmart(r(i-1),delta)**3)/3d0*
     *	(cos(theta(j-1))-cos(theta(j)))*4d0*pi
C- multiplied by 2 to account for 'southern' hemisphere
      	CM=CM+rho*dV
      	if(rho.gt.zero)then
      		V=V+dV
      	endif
      	CJ=CJ+rho*vphi*powrsmart(ra,delta)*sin(thetaa)*dV
      	CI=CI+rho*(powrsmart(ra,delta)*sin(thetaa))**2*dV
      P=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
      Pip=Pip+P*dV
      	T=T+rho*vphi**2/2d0*dV
      	W=W+rho*phia/2d0*dV
      	Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
      if((i.eq.IA).and.(j.eq.JMax))then
      omegaA=sqrt(max(zero,c1/(powrsmart(ra,delta)*sin(thetaa))))
      endif
      enddo
      enddo
      write(*,'(a,1p5g12.3)')' ra,thetaa,psi,c1,vphi=',ra,thetaa,psi,
     *	c1,vphi
C-        pause;
      return
      end
C---- '<--Leaving  Node %_integr:'
